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Abstract. We are developing a consistent method to combine tree-level event generators for hadron collision 
interactions with those including one additional QCD radiation from the initial-state partons, based on the 
limited leading-log (LLL) subtraction method, aiming at an application to NLO event generators. In this 
method, a boundary between non-radiative and radiative processes necessarily appears at the factorization 
scale (fip)- The radiation effects are simulated using a parton shower (PS) in non-radiative processes. It 
is therefore crucial in our method to apply a PS which well reproduces the radiation activities evaluated 
from the matrix-element (ME) calculations for radiative processes. The PS activity depends on the applied 
kinematics model. In this paper we introduce two models for our simple initial-state leading-log PS: a model 
similar to the "old" PYTHIA-PS and a pT-prefixed model motivated by ME calculations. PS simulations 
employing these models are tested using W^-boson production at LHC as an example. Both simulations 
show a smooth matching to the LLL-subtracted W + 1 jet simulation in the pr distribution of W bosons, 
and the summed pr spectra are stable against a variation of ^f, despite that the pr-prefixed PS results 
in an apparently harder pr spectrum. 
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1 Introduction 

Next-to-leading order (NLO) corrections in Quantum Chromodynamics (QCD) are required by experiments for precise 
measurements of known processes and reliable estimates of the background for expected heavy particle productions at 
high-energy hadron collisions, such as those at Fermilab Tevatron and the forthcoming CERN LHC. The corrections 
are desired to be implemented in Monte-Carlo event generators because they may depend on the topology of the events 
experimentally selected. 

The implementation of NLO corrections in event generators, in general, has a problem in regularizing the coUinear 
and soft divergences; i.e., the cancellation of divergences in radiative processes and non-radiative processes. Phase- 
space slicing methods that are usually applied in electroweak corrections are problematic for QCD corrections because 
of the large coupling. They frequently lead to large negative cross sections in non-radiative processes and/or visible 
boundaries in the phase space of radiative processes. These properties are undesirable for use in experimental analyses. 

For estimating interactions at high-energy hadron collisions, it is also necessary to employ an appropriate parton 
distribution function (PDF) and fragmentation functions (FF), depending on the energy scale of the perturbatively 
evaluated hard interaction. They resum large logarithmic collinear corrections to all orders to improve the convergence 
of the perturbation. These collinear corrections are usually simulated in event generators in the form of a parton 
shower (PS), because experiments require an exclusive simulation of beam collisions, and the corrections produce an 
experimentally visible transverse boost of the hard interaction and additional hadron activities. The corrections by 
PDF, FF and PS have an apparent overlap with real emissions in NLO corrections based on matrix-element (ME) 
calculations. We have to avoid double counts in order to construct consistent NLO event generators. 

Subtraction methods naturally solve these problems, in which the divergent terms are subtracted from radiative 
processes and their integrated contributions are added to non-radiative processes. The emergence of negative cross 
sections is limited and all results become finite, since the divergences are canceled within non-radiative processes. The 
application of PDF/FF and/or PS automatically restores the corrections, if the subtracted terms are identical to their 
leading contributions. The concept is simple, while implementation in event generators is not easy because the matching 
between the subtraction and the actual PS implementation is not trivial. After many trials, a practical solution [3] 
has been presented based on detailed knowledge on the initial-state PS of HERWIG [4|, and has been applied to 
several color-singlet production and heavy-flavor production processes where it is enough to consider radiations from 
initial-state partons. Recently, a new idea has been proposed to make the cross sections totally positive [5]. 

We are developing NLO event generators based on a different concept of subtraction [6]. Our goal is to develop 
an automated generation system of NLO event generators for hadron collision interactions based on GRACE [7] . We 
adopt a leading-log (LL) subtraction method for the matching between real emissions and PDF/PS, where a naive 
LL approximation to the parton radiation, proportional to is subtracted from radiative processes. Though the 

implementation is still limited to the initial-state radiation, extension to the final state will be easy since the concept 
is quite simple. 

In our method, radiative processes simulate the hard LL part, harder than the factorization scale {(J-f), and the 
non-LL part of the radiation. The softer LL part is simulated by the PS applied to non-radiative processes. Therefore, 
it is crucial in this method to employ a PS that well matches with the ME calculation of the radiation. In order to make 
the discussions transparent, we have developed an initial-state PS that strictly reproduces theoretical arguments in the 
LL approximation. A primary feasibility of this method has been demonstrated in a previous paper 1)] . The theoretical 
arguments are, however, performed only at the collinear limit. We need to introduce an appropriate kinematics model 
in order to construct a PS that is usable in practical event generators. The model applied in the demonstration was 
very primitive. 

In the present study we introduce two practical models of the PS kinematics and examine their matching properties, 
using W^-boson production at the LHC condition as an example. Since this process has an apparent fixed energy scale 
of the VT-boson mass (mw), possible mismatches would become visible in physical distributions. The matching can 
be tested without introducing loop corrections. The tests are carried out in the present study by using simulations 
employing tree-level ME calculations of the inclusive W (W + jet) and -I- 1 jet production processes, where "jet" 
denotes a light quark or gluon in the final state. We examine only the internal consistency of our method. Comparison 
with experimental data will be done in a separate paper, where additional lower effects, such as the intrinsic px 
and the underlying events, must be taken into account. 

We need to include loop corrections, together with cancelling soft and collinear corrections, in order to construct 
complete NLO event generators. Inclusion of these corrections will result in a modification to the normalization of 
non-radiative processes. Accordingly, it will produce a certain mismatch at the boundary {(J-f), even with a perfect 
matching method at the tree level. We will need to apply an appropriate modification to radiative processes, if we 
want to restore a smooth matching. By the way, this is merely a technical issue since this mismatch is a quantity 
at the level of NNLO (next-next-to-leading order), beyond the scope of NLO corrections. In addition we have a plan 
to construct a full NLL (next-to-leading log) parton-shower program [TU], in order to achieve a rigorous theoretical 
matching at the NLO level. The subtraction will have to be changed accordingly if it becomes available. 
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This paper is organized as follows: the bases of our LL-PS and the LL subtraction are given in Section 2. The 
description includes some details of the LL subtraction for the sample process. The kinematics models are introduced 
and tested in Sections 3 and 4. Finally, the conclusion is presented in Section 5. 



2 Parton Shower and LL subtraction 

We use an x-deterministic forward evolution [6j for the initial-state PS. This technique is a solution to overcome the 
low-efHciency problem in the forward evolution. Starting from a small Q^, the PS strictly reproduces the QCD evolution 
implemented in PDFs. The final Bjorken's x is thus given by the PS. The PS energy scale (/ips), the maximum energy 
scale in PS, is therefore necessarily identical to the factorization scale (hf), the energy scale of PDF. This means that 
the PS radiation is limited hy fip. The limitation is taken into account in the subtraction as well. Hereafter, we call 
this technique the limited leading-log (LLL) subtraction. 

In our PS we refer to a PDF at a low energy scale {Qq — 4.6 GeV) and produce branches in increasing order of 
until reaching = fij^. The of each branch is determined according to the Sudakov form factor, expressed as 



5(0?,Q^)=exp 
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where P{z) is the Altarelli-Parisi splitting function at the leading order; 
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for the branches q qg, g —>■ gg and g —>■ qq, respectively, with Cp = 4/3, Nq = 3 and Tr = nf/2 = 5/2. The 
splitting function P{z) is summed over possible branches. We use CTEQ5L [11] for the reference PDF in the present 
study. Accordingly we use the first-order strong coupling, 
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with /3o = 11 — 2n//3 and A — 0.146 GeV for nj ~ 5. The branching mode and the parameter z are determined in 
proportion to P{z) of relevant branches. An unphysical parameter, e, is introduced in order to cutoff the divergence 
of P{z). We set e = 10^^ as the default. Physical properties do not depend on this choice if e is set to be sufficiently 
small. 

PS is a Monte-Carlo solution of the QCD evolution based on the factorization theory. In this theory, the evolution 
is considered at the collinear limit in an infinite-momentum frame. The transverse behavior of the radiation is given 
only at the first-order approximation. As a result, the predicted kinematics do not strictly conserve the energy and 
momentum. We therefore need to introduce a certain model of the branching kinematics that gives a correspondence 
of the above branching parameters, and z, to kinematical variables in a finite-momentum frame, to construct a 
practical PS strictly conserving the energy and momentum. The purpose of this paper is to introduce such models 
and test their feasibility in our subtraction method. The use of our naive PS program allows us to clearly separate 
this model-dependent part from theoretically well-defined parts. 

Though the LL subtraction is not a main subject of this paper, we here describe some details for the present 
example process, -I- 1 jet production, since it is used for the tests of PS models in the following sections. At the 
parton level, the -I- 1 jet production process is composed of three incoherent subprocesses: 



qq 

g-q' 
qg 



>gW, 

>qW, 
q'W. 



(6) 

(7) 
(8) 



We ignore CKM nondiagonal interactions in the present study. Therefore, (17,9') denotes {d,u), {u,d), (s,c) or (c, s). 

The W bosons are assumed to decay to a pair of an electron and a neutrino. The matrix elements are evaluated 
including this decay at the tree level, assuming a total decay width of 2.12 GeV. The calculations are performed within 
the framework of BASES/SPRING ^ employing the matrix-element codes generated by GRACE [7]. The studies in 



4 



S. Odaka and Y. Kurihara: Initial-state parton shower kinematics for NLO event generators 




20 



40 



60 



-20 

100 



20 40 



60 



80 100 



p^(W) (GeV/c) 



Pt(W) (GeV/c) 



Fig. 1. -I- 1 jet cross section at LHC as a function of the VF-boson pr- (a) the sum of three incoherent subprocesses, and 
(b) the residual of the LL subtraction. A pt cut of 1 GeV/c is applied, while PS is yet to be applied in these results. The result 
from exact matrix elements at the tree level (filled circles) is compared with that from the LL approximation (histogram) in 
(a). The result is shown for each subprocess, qq' — > gW (solid), gq — > qW (dashed), and qg —* q'W (dotted), separately in (b), 
together with the sum shown with filled circles. 



this paper are carried out using cross-section results from BASES. The same method is applied for simulating the 
inclusive W {W + jet) production, too. 

Subtraction is done at the matrix- element level as [6] 

|A^LLsub(s;Aifl)P = - \Mw{.sw)? hi.{x,t] ^jlr), (9) 

where Aiw+j is the matrix element for a given -I- 1 jet event with a squared center-of-mass (cm) energy of s, and 
Aiw '^s that for the W production sub-system in this event having a squared invariant mass of s^. The rotation of 
the sub-system is also taken into account in the calculation, though it is irrelevant to the following studies since the 
angular properties are always integrated. The LL radiation factor /ll can be written as 

/LL(x,t,Mfl)-^^ (10) 

where P{x) is an Altarelli-Parisi splitting function relevant to this event. It is essential to use an identical definition 
of as in the -|- 1 jet matrix element M-w+j and in the LL factor /ll- We use the definition of Eq. (O with the 
energy scale set to a fixed value (^|.) in the present study. 

The parameter x is the ratio of the squared cm energies {x — sw / s) , and t is the squared momentum transfer of 
the parton branch in the picture of PS. There are two possibilities in the definition of t for the subprocess ([5]). The 
gluon in the final state can be radiated from q and q' in the initial state. We have to subtract both contributions. Thus, 
t is defined using the four-momenta of the partons as = l/|iq| -I- l/|iq'| with tq = {pq —Pg)^ and tqr — {pqi — pgY . 
On the other hand, there is no ambiguity for subprocesses ([7]) and ([H]), since g ^ qq is the unique possible branch. 
The parameter t is defined as t = {pg — pq)'^ for ([7]) and t = {pg — pq/)^ for dH). 

The subtraction result at the cross-section level is shown in Fig. [T] as a function of the transverse momentum of 
the W boson, where PS is yet to be applied. The matrix elements are converted to the cross sections assuming the 
LHC condition (proton-proton collisions at the cm energy of 14 TeV) and using CTEQ5L for PDF. The renormalization 
scale fin and the factorization scale hf are both fixed to the M^-boson mass {mw = 80.2 GeV). A px cut is necessary 
to apply since the subtraction is done numerically. We set the minimum pT to 1 GeV/c. These conditions are valid 
throughout the present study, except for the choice of fip- 
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Fig. 2. Factorization scale (/^f) dependence of the total cross section for W production at LHC. An apparent fip dependence 
in the inclusive W production cross section (open circles) is greatly reduced if the W + 1 jet cross section with the LLL (limited 
LL) subtraction is added (filled circles). Filled triangles illustrate the results from a simulation where the PS employing the 
PT-prefixed kinematics model is applied. The cross section result is affected by the PS veto in the VK + 1 jet simulation. The 
PDF is slightly changed since the QCD evolution is evaluated by PS in our simulation. 



Figure IHa) shows the cross section based on the exact + 1 jet matrix element (filled circles) and that based on 
the LL approximation (histogram); i.e., they correspond to the first and the second term in the right-hand side of Eq. 
([9]), respectively. The result from the subtracted matrix element is plotted in Fig.[TJb) with filled circles. The result is 
also shown for each subprocess separately. The residual is always negative for subprocess ([6]) and positive for ([7|) and 
([5]). The sum (filled circles) is positive since the contribution of the subprocess ^ is large in proton-proton coUisions. 
We can see that the residual is going to vanish as px goes to zero in all subprocesses, and the small pr cut does not 
significantly affect the final results. 

In the LLL (limited LL) subtraction, the subtraction is applied only when \tq\ < jj^p or \tqi \ < in the subprocess 
([U, and when \t\ < fj,"^ in the subprocesses ([7]) and ([5]). Here we note that both the -I- jet and the subtracted 
W -|- 1 jet cross sections depend on iip. The LL contributions of the parton radiation below the energy scale oi fip 
are resummed by PDF and PS in the former, while those above are explicitly evaluated at the first order of as in 
the latter. The non-LL contributions are all remaining in the latter. Thus, we can expect that the dependence on fj,p 
cancels in the sum of these two cross sections, at least at the leading order of ag. 

The summed total cross section is plotted in Fig. [2] with filled circles as a function of ^p- They should be compared 
with the inclusive W production cross section shown with open circles. We can see that the fj,F dependence is greatly 
reduced in the summed cross section, as we expected. The variation is at most 10% within the range of 1/2 < iip/mw < 
2. This is a remarkable feature of the LLL subtraction. The stability will become better if an appropriate resummation 
(Sudakov suppression) is applied to the 14^-1-1 jet process. It will suppress a small enhancement at very small fxp 
values in Fig.[2l as demonstrated in a previous paper [9]. 



3 PYTHIA kinematics 



We first test the kinematics model of the "old" PYTHIA-PS for the initial state [13] because this PS is constructed 
on nearly the same theoretical basis as ours. The starting assumption of this model is that the of the evolution is 
identical to the virtuality of the evolving parton. 



Pcvol — Q^' 



(11) 
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Fig. 3. Definition of kinematical variables in each parton branch. The calculation is done in the cm frame of the evolving 
parton and the target parton. 



and z gives the ratio between the squared cm energies after and before a branch, 



The latter ensures the relation Shard = X1X2S between the squared cm energies of the hard interaction and the beam 
collision, with x given by the product of all z values in each beam. 

The transverse momentum (pr) of each branch is calculated from energy-momentum conservation conditions |13[ 
[5]. The definition of kinematical variables is illustrated in Fig. [31 Assumption pT|) leads to = —Q^, while is 
given by the previous branch. The definition of z requires a target parton. The definition is not trivial since two 
incoming partons are both evolving. We execute the branches in the increasing order of equally taking both sides 
into consideration. The remaining parton on the opposite side is chosen as the target. Therefore, the target parton is 
also virtual. Going into further details, the calculation sometimes gives a negative p^. In such cases we set px to zero 
and adjust the to give this solution. These calculations are done in the cm frame of the evolving parton and the 
target parton as shown in Fig. [3l The momenta are rotated and boosted back to the cm frame of the starting partons 
every time after the branching kinematics are determined. After completing all branches, the system is boosted to 
the laboratory (beam-collision) frame, and the hard interaction part is attached according to the momenta of finally 
remaining two partons. 

Filled circles in Fig. HJa) show the simulated pr distribution of W^-bosons in the inclusive W production process 
for the LHC condition with /i_F = mw- The corresponding prediction from the "old" PYTHIA-PS is shown with a 
histogram, where PYTHIA 6.403 [E] is used with MSTP (68) = 0, MSTP(71) = 0, MSTP(81) = and MSTP(lll) = 0; 
i.e., the ME-correction, final-state PS, multiple-interaction and hadronization are turned off. The PDF choice is the 
default; namely, CTEQ5L is used, too. We can see a good agreement between these two simulations. Any difference in a 
very low px region (< 10 GeV/c) should be ignored because our Qo is relatively large and no intrinsic is considered. 
In the figure, PYTHIA predictions with the ME correction (MSTP(68) = 3) [M] using the "old" PS (MSTP(81) = 0) 
and the "new" PS (MSTP(81) = 20) are also illustrated for a reference. 

The inclusive W simulation result is also shown with a dotted histogram in Fig. [5] together with the result of 
the LLL subtracted ly -I- 1 jet simulation (dot-dashed histogram) described in the previous section, with the same 
energy-scale choice oi fip — mw- The PS is also applied to the -I- 1 jet simulation with a PS veto, where the PS 
is retried if the maximum pT of the PS branches exceeds the pr of the -I- 1 jet ME. The sum of these results is 
plotted with filled circles. We can see a smooth transition between the two simulations. As described in the previous 
section, the summed cross section shows good stability against the variation of ^p- A similar stability should also be 
observed in this pr spectrum if the applied PS kinematics shows a good matching to the ME calculation. 

The open circles and open squares in Fig.[5]show the summed results for other energy-scale choices of — 0.5mn/ 
and jip = l.5mw, respectively. We can see good stability. Though it is not shown in the figure, the result for 

— 2.0rn\Y shows a visible discrepancy from these results. This is reasonable since the PS is constructed on the basis 
of a theory and a model relevant to soft radiations. The underlying approximations may become inappropriate for 
hard radiations where the hardness {e.g., pr) is comparable to or larger than the typical energy scale of the considered 
hard interaction, mw in the present case. The contribution of such hard radiations may become significant if we set 
/x_F to be larger than mw- The energy scale fir should not be set very large in our method. 

Although the resultant simulation is self-consistent, the summed pr spectrum seems to be rather soft in low to 
medium regions (< 50 GeV/c) compared to the ME-corrected simulations of PYTHIA. This is because the applied 
PS is soft, as we can see in Fig. IH^a). Though this may be reasonable, we try to construct another kinematics model 
giving a harder spectrum in the next section. 

It may be worth noting here that the two ME-corrected PYTHIA simulations disagree with each other, as we can 
see in Figs. |4] and [3 The simulation with the "new" PS is harder and in good agreement with our simulation in the 



z = 



s'/s. 



(12) 
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Fig. 4. Simulated transverse momentum (pt) distribution of W bosons in the inclusive W production process at the LHC 
condition. Filled circles show the result of our simulation: (a) with the "old" PYTHIA-PS kinematics model and (b) with 
the pT-prefixed model. Solid histograms show the result from "old" PYTHIA without the ME correction to the initial-state 
radiation. The "old" and "new" PYTHIA simulations with the ME correction are illustrated with dashed and dotted histograms, 
respectively, for reference. Open triangles in (a) show the result of a pr-prefix model with a non-standard assumption of 
Py = (1 — z)^Q'^, and open circles in (b) the result of the "new" PYTHIA simulation without the ME correction. 



high-pT' region. Since the high-pr spectrum is predominantly determined by the unambiguous -I- 1 jet ME in our 
simulation, the "new" PS simulation must be reasonable. The softness of the "old" PS simulation that we see here 
may also be related to its basic property, which is discussed in the next section. 



4 pr-prefixed kinematics 

As we have shown in a previous report [9], the kinematics of the "old" PYTHIA-PS leads to the relation 

p| = (1-z)2q2 (13) 

in each branch at the soft Umit. As we can see in Eq. ([1]), the Sudakov form factor becomes smaller if we choose a 
smaller value for the unphysical parameter, e. This means that the number of parton branches increases. Thus, if we 
make e sufficiently small, the step of the PS evolution becomes nearly continuous compared to the energy scale of 
the hard interaction. However, the increased branches are predominantly very soft (z ~ 1) and do not affect physical 
properties, such as the pT of W bosons and visible jets. The above relation can be obtained at this continuous limit; 
i.e., p\-^ p^ — — using the notation in Fig. 

On the other hand, ordinary discussions concerning the branching kinematics give us a different relation, 

pI = {1- z)Q\ (14) 

This relation can be obtained by assuming that the incoming partons and the radiation are all massless; i.e. 
p1 = 0. Apparently Eq. (|13p gives smaller pt values than Eq. (|14p since z < 1. In ME calculations, all initial- and 
final-state partons are necessarily on-shell; namely, they are nearly massless. Therefore, if we want to achieve a good 
matching between the PS radiation and the ME radiation, the relation (fT4|) must be better than (|13p . We test such a 
model leading to the relation (fT4|) in this section. 

Although there may be some sophisticated models, here we adopt a straightforward way to realize this relation. We 
prefix the px of each branch according to Eq. (|14p as a starting assumption for calculating the branching kinematics. We 
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Fig. 5. Tree-level summation of the VK -I- jet simulation and the LLL-subtracted -I- 1 jet simulation for the LHC condition 
using the PS with the "old" PYTHIA-PS kinematics model. The W -(- jet result for the energy-scale choice of /xf = rnw is 
shown with a dotted histogram, the -1- 1 jet with a dot-dashed histogram, and the sum is plotted with filled circles. The 
summed results for fip = 0.5mw and /if ~ 1.5mw are plotted with open circles and open squares, respectively, in order to 
show the stability against the /iF variation. The dashed and solid histograms show the ME-corrected PYTHIA predictions for 
the W production using the "old" and "new" PS models, respectively. 



keep the definition of z, Eq. (|12[) . since otherwise we would need to apply undesirable energy/momentum corrections 
in the connection to ME. As a result, we have to abandon relation pip . The of the evolving partons is derived 
from energy-momentum conservation conditions. This must not be troublesome, since the virtuality is not a direct 
observable and the of the evolution is not a physically well-defined quantity, as we have already discussed. In the 
actual implementation, it sometimes happens that the prefixed pt exceeds the kinematically allowed maximum. We 
decrease the px value to the allowed maximum in such cases. 

The filled circles in Fig. HJ^b) show the result of a simulation with this new PS kinematics. The evolution to give 
and z is the same as in the previous simulation. Apparently it gives a harder px spectrum of W bosons. It should 
be noted here that relations and PH) are obtained at the soft limit {Q^ s). It is not trivial that the difference 
between these two relations is directly reflected in the visible px spectrum of W bosons. For a confirmation, we have 
tried another simulation with the pr-prefixed PS kinematics model, in which Eq. p^ is assumed instead of Eq. P4|) . 
The result is shown by the open triangles in Fig. [DJa). It is in good agreement with the previous result, confirming 
that the "old" PYTHIA-PS kinematics actually results in relation p3|) and the difference in the PS property at the 
soft limit is reflected in the visible px spectrum of W bosons. 

The matching test is retried using the new pT-prefixed PS in Fig. [6] The simulation conditions are the same as 
those giving the results in Fig. [5l except for the PS kinematics. We can see a smooth transition between the Vl^ -|- 
jet and 14^-1-1 jet simulations again. The spectrum in low to medium px regions has become harder and come closer 
to the PYTHIA simulations with the ME correction. The hf dependence is small, but has become larger as we can 
see in the figure. This must be simply due to the fact that this new PS generates harder radiations. Since the relations 
P3p and p4p are obtained at the limit where /s can be ignored, the approximation may become worse for hard 
radiations, as we have discussed in the previous section. In fact we can see a small bump structure in the fip /mw = 1 
result just below the boundary, though such a structure is not clear in the ^p/mw — 0.5 result. It is, however, 
possible to evaluate the internal ambiguity by comparing the predictions with different pip values; for instance, those 
with /ii? = my/ and = Q.bmw- In general, a Sudakov suppression has to be applied to the -I- 1 jet simulation 
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Fig. 6. Tree-level summation of the W -I- jet simulation and the LLL-subtracted W + 1 jet simulation for the LHC condition 
using the pT-prefixed PS kinematics model. The -)- jet result for the energy-scale choice of = rnw is shown with a dotted 
histogram, the IV -I- 1 jet with a dot-dashed histogram, and the sum with filled circles. The summed results for fip = 0.5mw 
and fiF = 1.5mw are plotted with open circles and open squares, respectively. The dashed and solid histograms show the 
ME-corrected PYTIfIA predictions for the W production using the "old" and "new" PS models, respectively. 



in order to take higher-order effects into account, if we choose an energy scale apart from the typical energy scale of 
the considered hard interaction. The suppression is however a few percent and can be ignored for — Q.^mw- 

It should be noted here again that the Bjorken's x for the hard interaction is given by PS in our method. The PS 
veto, therefore, alters the PDF for the IV -|- 1 jet simulation. The cross section prediction slightly changes as a result. 
The filled triangles in Fig. [2] show the results of a simulation with the pT-prefixed PS assuming the relation At 
present we have no idea about the question whether this alternation is reasonable or not. 

In the course of the present study, we have found that the "new" py-ordered PS of PYTHIA shows a behavior 
quite similar to our pr-prefixed PS. The result of the "new" PYTHIA-PS simulation without the ME correction is 
over-plotted in Fig. [4l[b) with open circles. We can see that the pt spectrum is almost identical to our new result. 



5 Conclusion 

In the limited leading-log (LLL) subtraction method for NLO event generators, the parton shower (PS) to be applied to 
non-radiative processes must well reproduce the leading-log (LL) contributions in the matrix-element (ME) evaluation 
of radiative processes. The PS activity depends on the applied model for branching kinematics. We have introduced 
two kinematics models to be applied to our simple initial-state leading-log PS: a model similar to the "old" PYTHIA- 
PS and a p-r-prcfixed model with a pt definition motivated by ME calculations. The former shows a behavior similar 
to the "old" PYTHIA-PS, and the latter results in harder transverse activities, as we expected. We also found that 
the "new" PYTHIA-PS is similar to the latter. 

PS simulations employing these two models have been tested using the W^-boson production at LHC as an example. 
Both simulations show a smooth matching to the LLL subtracted -I- 1 jet simulation in the pT distribution of W 
bosons, though the pT-prefixed PS results in an apparently harder pT spectrum than the PS with the "old" PYTHIA 
kinematics. The summed pr spectra are stable against the variation of the factorization-scale choice in both simulations. 
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This is a remarkable feature of the LLL subtraction method. The remaining instabihty can be used as a measure of 
the internal ambiguity of this method. 

The present study is not enough to decide which model is better for our application. The decision must be postponed 
until further studies can be conducted employing a comparison with experimental data, such as the .Z-boson production 
data at Tevatron. Lower effects that are absent from the present simulations has to be taken into account there. 

The method tested in this paper is not dedicated to M^-boson production, but can be apphed to other processes. Of 
course, my/ has to be replaced with their typical energy scales in such applications. At present, our method is limited 
to applications to the initial-state radiation. It will, however, be easy to extend it to the final state since the concept is 
very simple. Once extended, it will allow us to construct NLO event generators for those processes including jet(s) in 
the final state, which the existing NLO event generator [3J has not yet supported. Even if loop corrections necessary 
for NLO calculations are absent, our subtraction method will give us a consistent way to simulate multi-jet production 
processes at the tree level. This will be an alternative to the CKKW method [15] if it can be applied recursively. 
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